source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_result.R')
library(hmod)
library(magrittr)
library(data.table)
library(lattice)
library(scatterplot3d)
library(microbenchmark)
library(ggplot2)
library(grid)
library(gridExtra)
library(pracma)
library(here)
text.size = 10
dir = here()
# energy firm size regression
############################################
wd = gsub("Figures", "Empirical Data/GEM", dir)
setwd(wd)
energy = fread("energy_v_firm.csv")
e_predict = function(firm_size){
x = log(energy$energy)
y = log(energy$firm_mean)
r = lm(y ~ x)
a = exp(coef(r)[1])
b = coef(r)[2]
y = (firm_size/a)^(1/b)
return(y)
}
# model
######################################################
wd = "/home/blair/Desktop/within_firm/Supplementary Material/hierarchy_model/data"
setwd(wd)
compustat = fread("compustat.txt")
s = fread("s_empirical.txt")
g = fread("g_empirical.txt")
span = boot_span(s$V1, s$V2)
a = span[1]
b = span[2]
sigma = boot_sigma(g$V1)
# tune to compustat
base_compust = base_fit(a, b, compustat$V1)
fit = fit_model(a, b, base_compust, compustat$V1, compustat$V2, compustat$V3, 0.001)
fit = data.table(compustat, fit)
fit = fit[fit$ceo_error < 0.001]
#  low energy
##########################################################################################
# simulation
wd = paste(dir, "/video_test", sep = "")
setwd(wd)
n = 10^4
n_iterations = 100
firm_sample = rpld(10^7, 1, 3.5, 2.3*10^4)
growth_rate = rnorm(n, 1, 0.5)
hist(growth_rate)
library(truncnorm)
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.5)
hist(growth_rate)
firm_new = firm_sample*growth_rate
firm_new = round(firm_sample*growth_rate)
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
low_mean_firm = mean(firm_new)
low_energy = e_predict(low_mean_firm) %>% round()
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
mean(firm_sample)
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.2)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.2)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.2)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.2)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.2)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.2)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.2)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.2)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.2)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.2)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.2)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.2)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.2)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.2)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
mean(firm_sample)
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.3)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.3)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.3)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.3)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.3)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.3)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.3)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.3)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
library(snow)
library(doSNOW)
library(foreach)
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
stopCluster(cl)
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.3)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
low_mean_firm = mean(firm_new)
low_energy = e_predict(low_mean_firm) %>% round()
firm_sim = sample(firm_sample, 10^6)
base_sim = base_fit(a, b, firm_sim)
r_sim = r_sim_2(fit$V1, fit$r, firm_sim)
base_pay_sim = bp_sim(fit$base_pay, 10^6)
mod_low = model(a, b, base_sim, firm_sim, base_pay_sim, r_sim, sigma, power = T) %>% data.table()
low_power_gini = fastgini(mod_low$power) %>% round(.,2)
low_energy_label = paste(  "Energy Use per Capita ~ ",
low_energy, "GJ per person",
"\nGini Index of Hierarchical Power Concentration =", low_power_gini)
# firm grid
###########################################################################################
firm_sim = sample(firm_sample, n)
base_sim = base_fit(a, b, firm_sim)
r_sim = r_sim_2(fit$V1, fit$r, firm_sim)
base_pay_sim = bp_sim(fit$base_pay, n)
coord = firm_grid(base_sim, 20)
grid = grid_plot(a, b, sigma, coord, base_sim, base_pay_sim, r_sim)
colnames(grid) <- c("Firm_ID", "x_", "y_", "z_", "pay")
# camera parameters
e = c(-4, -4, -4)
c1 = c(-300, -300, 400)
theta = c(pi/2, 0, -pi/2.2)
# 2D projection
r = ((grid[,2]-c1[1])^2 + (grid[,3] - c1[2])^2)^0.5
grid = grid[order(-r),]
p.bottom = cbind(grid[, 2:3], z_ =  exp( 0.4*(grid[,4] - 1) ) - 1 )
bottom = project(p.bottom, c1, theta, e)
colnames(bottom) <- c("x", "y")
grid[,4] =  exp( 0.4*grid[,4])
top = project(grid[, 2:4], c1, theta, e)
colnames(top) <- c("x", "y")
people.coord = data.table(bottom = bottom, top = top, pay =  grid[,4])
x_max = max(coord[,3])
y_max = max(coord[,5])
border = rbind(
c(-1, -1, 0),
c(x_max + 1, 0, 0),
c(x_max + 1, y_max + 1, 0),
c(-1, y_max + 1, 1)
)
border.coord = project(border, c1, theta, e) %>% data.table()
g_low = ggplot() +
geom_polygon(data = border.coord, aes(x = V1, y = V2), fill = NA, col = "grey60", size = 0.1) +
geom_segment(data = people.coord,
aes(x = bottom.x, y = bottom.y, xend  = top.x, yend = top.y, col = pay), size = 0.2) +
scale_color_gradientn("Relative Income\n(median = 1)\n",    colours = rainbow(4)) +
coord_cartesian(xlim = c(min(people.coord$bottom.x), max(people.coord$bottom.x))*0.95, ylim = c(-1, 2.5)) +
theme_bw() +
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
legend.key.width = unit(0.3, "cm"),
legend.text = element_text(size = rel(0.8)),
legend.title = element_text(size = rel(0.8), vjust = 0),
text=element_text(size= text.size, family="Times"),
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank()
) +
annotate("text", label = "Subsistence Society",
x = mean(range(people.coord$bottom.x)) , y = 2.5, family = "Times", size = 4.5, fontface = 2 ) +
annotate("text", label = low_energy_label,
x = mean(range(people.coord$bottom.x)) , y = 2, family = "Times", size = 3.5 )
stopCluster(cl)
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
mean_firm = mean(firm_new)
energy = e_predict(mean_firm) %>% round()
firm_sim = sample(firm_sample, 10^6)
base_sim = base_fit(a, b, firm_sim)
r_sim = r_sim_2(fit$V1, fit$r, firm_sim)
base_pay_sim = bp_sim(fit$base_pay, 10^6)
mod_low = model(a, b, base_sim, firm_sim, base_pay_sim, r_sim, sigma, power = T) %>% data.table()
power_gini = fastgini(mod_low$power) %>% round(.,2)
energy_label = paste(  "Energy Use per Capita ~ ",
energy, "GJ per person",
"\nGini Index of Hierarchical Power Concentration =", power_gini)
# firm grid
###########################################################################################
firm_sim = sample(firm_sample, n)
base_sim = base_fit(a, b, firm_sim)
r_sim = r_sim_2(fit$V1, fit$r, firm_sim)
base_pay_sim = bp_sim(fit$base_pay, n)
coord = firm_grid(base_sim, 20)
grid = grid_plot(a, b, sigma, coord, base_sim, base_pay_sim, r_sim)
colnames(grid) <- c("Firm_ID", "x_", "y_", "z_", "pay")
# camera parameters
e = c(-4, -4, -4)
c1 = c(-300, -300, 400)
theta = c(pi/2, 0, -pi/2.2)
# 2D projection
r = ((grid[,2]-c1[1])^2 + (grid[,3] - c1[2])^2)^0.5
grid = grid[order(-r),]
p.bottom = cbind(grid[, 2:3], z_ =  exp( 0.4*(grid[,4] - 1) ) - 1 )
bottom = project(p.bottom, c1, theta, e)
colnames(bottom) <- c("x", "y")
grid[,4] =  exp( 0.4*grid[,4])
top = project(grid[, 2:4], c1, theta, e)
colnames(top) <- c("x", "y")
people.coord = data.table(bottom = bottom, top = top, pay =  grid[,4])
x_max = max(coord[,3])
y_max = max(coord[,5])
border = rbind(
c(-1, -1, 0),
c(x_max + 1, 0, 0),
c(x_max + 1, y_max + 1, 0),
c(-1, y_max + 1, 1)
)
border.coord = project(border, c1, theta, e) %>% data.table()
g_low = ggplot() +
geom_polygon(data = border.coord, aes(x = V1, y = V2), fill = NA, col = "grey60", size = 0.1) +
geom_segment(data = people.coord,
aes(x = bottom.x, y = bottom.y, xend  = top.x, yend = top.y, col = pay), size = 0.2) +
scale_color_gradientn("Relative Income\n(median = 1)\n",    colours = rainbow(4)) +
coord_cartesian(xlim = c(min(people.coord$bottom.x), max(people.coord$bottom.x))*0.95, ylim = c(-1, 2.5)) +
theme_bw() +
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
legend.key.width = unit(0.3, "cm"),
legend.text = element_text(size = rel(0.8)),
legend.title = element_text(size = rel(0.8), vjust = 0),
text=element_text(size= text.size, family="Times"),
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank()
) +
annotate("text", label = "Subsistence Society",
x = mean(range(people.coord$bottom.x)) , y = 2.5, family = "Times", size = 4.5, fontface = 2 ) +
annotate("text", label = energy_label,
x = mean(range(people.coord$bottom.x)) , y = 2, family = "Times", size = 3.5 )
# export
#############################################################################
gA = ggplotGrob(g_low)
png(sprintf("frame%05d.png", i) , width = 8, height = 4, units = 'in', res = 600)
grid.arrange(gA, nrow = 1)
dev.off()
stopCluster(cl)
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
i = ``
i = 1
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.1)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
mean_firm = mean(firm_new)
energy = e_predict(mean_firm) %>% round()
firm_sim = sample(firm_sample, 10^6)
base_sim = base_fit(a, b, firm_sim)
r_sim = r_sim_2(fit$V1, fit$r, firm_sim)
base_pay_sim = bp_sim(fit$base_pay, 10^6)
mod_low = model(a, b, base_sim, firm_sim, base_pay_sim, r_sim, sigma, power = T) %>% data.table()
mod = model(a, b, base_sim, firm_sim, base_pay_sim, r_sim, sigma, power = T) %>% data.table()
power_gini = fastgini(mod$power) %>% round(.,2)
energy_label = paste(  "Energy Use per Capita ~ ",
energy, "GJ per person",
"\nGini Index of Hierarchical Power Concentration =", power_gini)
firm_sim = sample(firm_sample, n)
base_sim = base_fit(a, b, firm_sim)
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.3)
firm_new = round(firm_sample*growth_rate)
# stochastic growth
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.3)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
mean_firm = mean(firm_new)
energy = e_predict(mean_firm) %>% round()
firm_sim = sample(firm_sample, 10^6)
base_sim = base_fit(a, b, firm_sim)
r_sim = r_sim_2(fit$V1, fit$r, firm_sim)
base_pay_sim = bp_sim(fit$base_pay, 10^6)
mod = model(a, b, base_sim, firm_sim, base_pay_sim, r_sim, sigma, power = T) %>% data.table()
power_gini = fastgini(mod$power) %>% round(.,2)
energy_label = paste(  "Energy Use per Capita ~ ",
energy, "GJ per person",
"\nGini Index of Hierarchical Power Concentration =", power_gini)
energy = e_predict(mean_firm) %>% round()
e_predict(mean_firm)
mean_firm = mean(firm_new)
energy = fread("energy_v_firm.csv")
# energy firm size regression
############################################
wd = gsub("Figures", "Empirical Data/GEM", dir)
setwd(wd)
energy = fread("energy_v_firm.csv")
stopCluster(cl)
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
stopCluster(cl)
i = 1
firm_sample = rpld(10^7, 1, 3.5, 2.3*10^4)
d
5
stopCluster(cl)
wd = "/home/blair/Desktop/within_firm/Supplementary Material/hierarchy_model/data"
setwd(wd)
compustat = fread("compustat.txt")
s = fread("s_empirical.txt")
g = fread("g_empirical.txt")
span = boot_span(s$V1, s$V2)
a = span[1]
b = span[2]
sigma = boot_sigma(g$V1)
# tune to compustat
base_compust = base_fit(a, b, compustat$V1)
fit = fit_model(a, b, base_compust, compustat$V1, compustat$V2, compustat$V3, 0.001)
fit = data.table(compustat, fit)
fit = fit[fit$ceo_error < 0.001]
# simulation
wd = paste(dir, "/video_test", sep = "")
setwd(wd)
n = 10^4
n_iterations = 100
firm_sample = rpld(n, 1, 3.5, 2.3*10^4)
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.3)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sample = firm_new
mean_firm = mean(firm_new)
energy_mod = e_predict(mean_firm) %>% round()
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
firm_sim = rpld(n, 1, 3.5, 2.3*10^4)
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.3)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
firm_sim = firm_new
mean_firm = mean(firm_new)
energy_mod = e_predict(mean_firm) %>% round()
base_sim = base_fit(a, b, firm_sim)
r_sim = r_sim_2(fit$V1, fit$r, firm_sim)
base_pay_sim = bp_sim(fit$base_pay, 10^6)
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
growth_rate = rtruncnorm(n, a = 0, b = Inf,  1, 0.3)
firm_new = round(firm_sample*growth_rate)
firm_new [ firm_new == 0] = 1
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
1.1^100
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_landscape_video.R')
